Stress Propagation in Two Dimensional Frictional Granular Matter 
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The stress profile and reorientation of grains, in response to a point force applied to a preloaded 
two dimensional granular system, are calculated in the context of a continuum theory that incor- 
porates the texture of the packing. When high friction prevents slip at the inter-grain contacts, an 
anisotropic packing propagates stress along two peaks which amalgamate into a single peak as the 
packing is disordered into a less anisotropic structure; this single peak may be wider or narrower 
than in the isotropic case, depending on the preparation of the packing. At lower frictions, an 
effective treatment of slipping contacts yields sharpened peaks, and ultimately a singular limit in 
which stress propagates along straight rays. Recent experiments, as well as aspects of hyperbolic 
models, are discussed in light of these results. 



Descriptions of static granular matter as an isotropic 
elasto-plastic continuum |jj do not incorporate effects of 
granularity essential to reproduce the variety of stress 
profiles identified experimentally 0-||]. By contrast, re- 
cent discrete models relying on (probabilistic) rules of 
force transmission among individual grains |^,|l^,|ll| re- 
quire, for a continuum limit in terms of stresses, an ad- 
ditional constraint relating the components of the stress 
tensor. This closure relation may be interpreted |TT[ as a 
local Janssen approximation |l^,[lj but is otherwise rather 
ad hoc. Furthermore, it leads to hyperbolic equations for 
stress propagation Jl3|,[l4|,[ll| , as opposed to the usual el- 
liptic equations of elastic theories. 

Reference offers a possible bridging approach in 
which a continuum theory is derived systematically from 
grain-grain interactions, taking into account the texture 
of the packing and the freedom of grains to reorient with 
respect to their surrounding medium. In the present Let- 
ter, we apply this general framework to calculate the 
propagation of the stress from a point force applied to 
a two dimensional granular system. After recalling the 
main elements of the theory and discussing its general 
properties in two dimensions, we calculate stress profiles 
in the presence of a friction large enough to prevent slip 
at the contacts. We obtain a variety of responses that 
depend on the type of packing and provide us with a 
natural interpretation of recent experiments. Finally, we 
discuss the case of easily sheared contacts, which mimic 
slipping contacts and induce sharper stress peaks. In the 
limiting situation, the hyperbolic closure relations used 
in the literature [ fH3|Ji!|Ji"l| ] emerge from our formulation 
and constrain the applied point force to propagate along 
straight rays. 

As shown in Ref. a strain e about a preloaded 

equilibrium state induces reorientation of the grains, ex- 
pressed by the linearized rotation matrix Sij 
generates incremental stresses 
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encode the texture of the granular packing via the aver- 
age number of inter-grain contacts fi(a)da within an an- 
gle da about a, the normal and shear moduli fen (a) and 
k±(a) of these contacts, and the average center-to-center 
distances D(a) traversing them. In two dimensions O 
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and set by the local vanishing of torque to 



depends on a single angle uj, defined by Q = 
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The stress then can be expressed in terms of the strain 
alone, for example in the matrix form 



M 
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where the (symmetric) compliance matrix M summa- 
rizes the textured stiffness of the medium through 
the first few Fourier moments of the quanti- 
ties fJ-Q{a) = k±(a) D 2 (a)[i{a) and fip(ct) — 
[fcii (a) — fcj_(a)] D 2 (a)ii(a) associated with the two fab- 
ric tensors. Defining the Fourier moments by 



2-k^q{o) — 2k + 4/ti cos(2a) + 4k 2 sin(2a) 
2TT/j,p(a) = 2R + 4Ri cos(2a) + 4k 2 sin(2ci!) 

+ I6K3 cos(4a) + I6R4 sin(4a) + 
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we write M explicitely, as 

K + + K 3 — K 3 k 2 + K 4 

M = I —K3 k - Kl + K 3 k 2 - K 4 J , (7) 

k 2 + k 4 k 2 - ki — R3 + k 5 

where k = k + R , k\ = K\ + R\, k 2 = (k 2 + R 2 )/2, 
ks = R3 — So/4 — k 2 /2kq, K4 = K4 + kik 2 /2kq, and 
«5 = kq/2 — (k\ + k 2 )/2kq. The fore-aft symmetry 
(a — > a + 7r) of the contact distribution forbids dipo- 
lar moments, and the anisotropic part of M is composed 
of the quadrupolar moments Ki, Ri, n 2 , R 2 and octupolar 
moments S3, R4. 

Before discussing stress profiles generated by Eq. (||), 
we comment on the properties of the compliance. The 
non-negativity of [aq and /ip restricts the allowed domain 
of the multipolar components; Rq, R\ and S3, for example, 
are bound by 2(ki/So) 2 — 1 < AR^/Rq < 1. An intuitive 
derivation of these (Cauchy-Schwartz) inequalities uses 
a representation of fiQ and \xp as linear combinations of 
delta functions with non-negative weights. The multipo- 
lar components of [iq and /ip then are themselves linear 
combinations of the multipolar components of the delta 
functions with the same non-negative weights, hence the 
former interpolate the latter. For a (fore-aft symmetric) 
pair of delta functions /i ao (a) = 5(a — ao) + 5(a — ao — ir), 
ko(Ha ) = 1, KiOO = cos(2a ), k 3 (n ao ) = cos(4a )/4, 
and the bound 2(ki/So) 2 — 1 < AR^/Rq is achieved. The 
choice fi(a) — fJ. ao { a ) corresponds to the rather patho- 
logical limit of a granular medium made up of an array 
of independent, parallel columns at an angle ao, and lies 
on the boundary of the allowed domain. This domain 
is included entirely in the domain of stability of M , in 
which stress resists any deformation (with a positive en- 
ergy cost). As long as friction is high enough to prevent 
slip at the contacts, the boundaries of the domain of sta- 
bility and the allowed domain meet only at the columnar 
systems just mentioned. Finally, the form of M in Eq. ([?]) 
and the independence of Hq, . . . , K5 ensure that, within 
the allowed domain, the granular material can take on 
the most general form of anisotropic elasticity. 

The equilibrium conditions dj <Jij — fi S(x)5(y), where 
fi is a point force applied at the origin, are supplemented 
by the closure equation d xx e yy + d yy e xx = 2d xy e xy . This 
relation specifies the constraint on a strain tensor that 
derives from an underlying two-component displacement 
vector field, and is re-expressed as 



{d yy , d x 
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to close the system of equations governing stresses. A 
classic and elegant solution makes use of the Airy func- 
tion <j), defined through its derivatives by Sijd KK <p — 
didj(f> = <Jij so as to satisfy the equilibrium conditions 
automatically. Thus one needs to solve only Eq. (pi) for 



respecting the boundary conditions at the origin and 
at infinity. Taking advantage of the properties of analytic 
functions on the plane, one can translate this fourth order 
differential equation for <j) into a fourth order algebraic 
equation Jl6[ . Its roots determine completely the angular 
modulation that factors the usual 1 jr dependence of two 
dimensional elastic stresses. In an orthotropic material 
(with two perpendicular axes of symmetry), the quar- 
tic polynomial in question reduces to a bi-quadratic one 
Jl7| , and henceforth we restrict ourselves to this case for 
the sake of simplicity. Specifically, we choose the x- and 
y-axes as axes of symmetry, so that k 2 — R 2 = K4 = 0. 

Following the program just outlined, we calculate the 
radial stress a rr and the reorientation angle oj analyti- 
cally and plot their angular dependence on Fig. 1 p8| . 
While observed modulated ( e.g., double-peaked) re- 
sponses motivate hyperbolic models Jl^JlJ,|lJ], it was em- 
phasized recent ly |20| that they arise also in anisotropic 
elastic systems fl7j , |l6 ,^l|. Accordingly, for a weakly dis- 
ordered square lattice oriented at 45° from the x- and 
y -axes (R3 < 0), we find that the point force propa- 
gates through the material as a shallow stress double- 
peak (middle bold curve in Fig. 1(a)). As the packing is 
further disordered, the two peaks merge into a single peak 
which ultimately coincides with the isotropic response 
(thin solid curve). Alternatively, if contacts close to the 
vertical are slightly favored (k\,Ri < 0) e.g., by gravity, 
the two peaks amalgamate into a single wide peak (top 
bold curve). If, however, the preparation of the system 
favors contacts close to the horizontal (kx,Ki > 0), the 
two peaks become more pronounced (bottom bold curve) 
and the stress propagates mostly along two directions. 
(We note that there is no trivial relation between the 
latter and possible preferred directions in the packing.) 
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1. Radial stress and reorientation angle as a func- 
tion of the polar angle, for different packings of incompress- 
ible grains (ki = 2R\ [ p^|Jl^ ]). (a) Bold curves: K3 = —1/5; 
/ti = —1/5 (top), Ri = (middle), Ri = 1/5 (bottom). Dot- 
ted curve: Ri — —7/10, ^3=0. The isotropic result (thin 
solid line) is added for comparison, (b) Ri — 1/5, R3 = —1/5. 
The crosses illustrate the angular dependence of the reorien- 
tation and its 1/r decay. 



The plots are symmetric because they correspond to 
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packings statistically symmetric about the vertical (di- 
rection of the force). If the material is rotated with re- 
spect to the force, one of the two peaks dominates the 
other; more generally, if the force is not applied along 
an axis of symmetry (as generically in a non-orthotropic 
medium), we expect a biased response reminiscent of 
arching phenomena |2^ ]. This is certainly the case in 
the left (or right) half of a pile constructed with a point 
source, where the packing is not symmetric about the ver- 
tical dJl^l- When K\ 7^ 0, the stress profile is due in part 
to the reorientation of the grains. Figure 1(b) exhibits 
the angular dependence of the reorientation angle ui cor- 
responding to the pronounced double-peak in Fig. 1(a); 
the crosses illustrate the form of the reorientation field 
and its decay away from the origin. While the reorienta- 
tion appears to relieve the shear stress along the vertical 
as expected intuitively, it shows a complicated structure 
and in particular extrema along non-trivial directions. 

In a remarkable measurement of the response of a 
vertical two dimensional granular system (preloaded by 
its own weight) to a vertical point force, the stress propa- 
gates along two peaks, in an ordered (anisotropic) piling, 
which merge into a single peak as the piling is disor- 
dered, in conformity with our picture. Furthermore, the 
width of the single peak increases linearly with depth, 
in accordance with an elastic theory. (The dependence 
of the widths on depth in the case of a double-peak is 
not given.) Another experiment || correlates the macro- 
scopic response to the microscopic structure by record- 
ing both the stress profile and the angular distribution 
of contacts /i(a). As expected fJ,(a) shows a higher- 
degree of anisotropy in a pile grown with a point source 
than in a pile grown with an extended source, and dif- 
ferent stress profiles are measured: a double-peak for a 
point source, a single peak for an extended source. In 
similar experiments on three dimensional piles ||[| the 
observed stress peak, while scaling linearly in depth, is 
narrower than expected from isotropic elasticity in an in- 
finite half space B. The authors interpret this accentu- 
ated response as possibly arising from the experimental 
boundary conditions |^] — a conjecture left unconfirmed 
by a more detailed study M. Our approach affords a 
complementary explanation: no matter how one disor- 
ders the pile during or after its growth, gravity sets a 
preferred direction; one then expects to have more near- 
vertical contacts than near-horizontal ones, and a result- 
ing narrower response (dotted curve in Fig. 1(a)). This 
interpretation is supported by measurements Q showing 
that a dense packing yields a wider peak than a loose 
packing. To go from dense to loose, one needs to "open 
holes" while maintaining stability; thus these holes may 
decrease the number of near-horizontal contacts relative 
to the number of near-vertical ones, but not the other 
way around. 

To complete our picture of stress propagation, we 
discuss the degenerate solutions that arise when slip 



occurs. As long as a large friction prevents slip, 
the ratio of shear to normal moduli of the contacts 
is fixed (for purely compressively preloaded contacts, 
k±/k\\ = (2 — 2v) I (2 — v), where v is the Poisson ratio 
of the constitutive material of the grains |l9|]). At smaller 
frictions, grains may slip along their contacts, effectively 
lowering the value of k± p3]. This renormalization of 
the shear modulus yields narrower stress peaks, and in 
the limit k±/k^ — > singular responses arise in ordered 
packings. For a simple two dimensional crystal (Bravais 
lattice) fj, Q (a), fip(a) oc X)"=i^( a - oti) + 5(a - a t - tt)], 
where steric constraints impose n < 3 24 1. The case 
n = 1 corresponds to the pathological columnar systems 
mentioned above, while n — 3 yields an isotropic elas- 
tic theory since the contact angles differ by multiples of 
7r/3. In the case n — 2, the stress peaks become in- 
finitely narrow and coincide with the lattice directions 
«i and c*2- Similar ray- like responses have been observed 
recently in three dimensional crystalline packings Q . In 
two dimensions, the situation is easiest to visualize in a 
square lattice at 45° (Fig. 2); a point force applied ver- 
tically is transmitted along two rows of grains without 
disturbing any other grains. We note that, in general, 
in the limit k_\_/k\\ — > only the fabric tensor P sur- 
vives in Eq. (|l]) and the system becomes equivalent to 
a collection of fibers highly resistant to compression and 
extension but easy to shear, as in much studied fiber- 
reinforced materials |25| . Ordering corresponds to align- 
ing the fibers along a set of given directions; in the case 
of only two directions there is at least one soft mode, 
associated with a singular response. In particular, the 
conformation of Fig. 2 is equivalent to two sets of fibers 
at right angle, equivalent in turn to an incompressible 
elastic medium with fibers along a single direction. 




FIG. 2. Illustration of a singular response in an ordered 
packing. Left: the force is transmitted through the two rows 
of black grains. Right: detail; over- and under-compressed 
contacts are in black, unperturbed ones in white. 



Singular solutions of a hyperbolic character occur only 
at the boundary of the domain of stability of the elas- 
tic theory. In the limit k±/k^ — > 0, this boundary and 
that of the allowed domain meet not only at the n = 1 
columnar systems but also at the n = 2 crystalline pack- 
ings. For a packing symmetric with respect to the x - 
and y-axes (ai = — «2 = cvq), the closure Eq. <M) reduces 
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to 



where c, 



(d xx 



C 9yy)(a xX C (Tyy) — 0, 



(9) 



(«o - «i)/(«o/4 — K3) — 1 = tan 2 a > 0. If 
the system is rotated by an angle arctan(t), the closure 
equation involves the more general linear combination 
(l - cgi 2 ) a xx - (c 2 , - t 2 ) cj y y - 2(1 + c 2 ,) ta xy . Since the 
limit A;j_/fc|| — > may be viewed also as that of infinitely 
hard grains (with a thin incompressible elastic coating), 
it is no surprise that we retrieve the form of the closure 
equations derived |26|,^7j for the marginal (isostatic) state 
of a granular packing in which the response may be cal- 
culated without knowledge of deformations. In integral 
form, our closure equation with Cq = 1 is identical to 
the one given in Ref. while for general Co and t it 
coincides with that of Ref. f2^J and the phenomenologi- 
cal "oriented stress linearity model" of Ref. |l4|] . In Ref. 
jnj, an integral form of Eq. @ is suggested by a dis- 
crete probabilistic model and the packing disorder is em- 
ulated by adding to c% a small increment that fluctuates 
randomly in space. This uncertainty in the local propa- 
gation directions results in a diffusive broadening of the 
stress rays, proportional to the square root of the depth. 
Within our approach, by contrast, as soon as disorder 
in the contact orientation is introduced, the response of 
a homogeneous system is governed by elastic (elliptic) 
equations on large scales, yielding a linear broadening 
of the stress peaks. These remarks help conciliate the 
square-root broadening observed Q| in a pile less than 
ten grains in height with the linear broadening J|-|7j on 
larger scales. The situation might be different, however, 
for heterogeneous media; in a polycrystalline packing, for 
example, one might expect a diffusive behavior even on 
macroscopic scales. 
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